!-----------------------------------
! A. Gazizov, LNGS, 18.06.11

module coeffnoosc
   implicit none

   private AprBD

! ee-pair production on CMB
   integer,parameter,private        &
       :: nee   = 201                 ! number of energy lines	   
   real(8),dimension(nee),private   &
       :: lgee, lgbetaee &	          ! lg(E/eV), lg(beta^{ee}/H_0)
        , dbetaee	  		          ! d\beta^{ee}/dlnE

! N + \gamma -> N' + X on CMB 
   integer,parameter,private        &
       :: necmb = 326                 ! number of energy lines
   real(8),dimension(necmb),private & 
       :: lgecmb                      ! log10(E/eV)

! The case of quick oscillation
! No neutron components are to be taken into account
   real(8), private:: tt   ! foo
   real(8), dimension(necmb), private   &
      :: Ppn                            &
       , Bpp, dBpp, Dpp, dDpp, d2Dpp    & 
       , Bpn, dBpn, Dpn, dDpn, d2Dpn 
   real(8), dimension(13), private:: C

   real(8), dimension(5):: Cour	 ! output of this module


!------------------------
         CONTAINS
!------------------------

!-----------------------------------------
   Subroutine ReadPCMB
      implicit none
	  integer:: j

! ------------ Arrays for CMB  ---------------
      open(1,file='Data/BBPee.dat',status='old')
	     read(1,*)
	     read(1,*)
         read(1,'(f8.4,f18.11,E20.11)')               &
            (lgee(j),lgbetaee(j),dbetaee(j), j=1,nee)
	  close(1)
!------------- 
	  open(1,file="Data/pwij.dat",status='old')
	     read(1,'(//(25E15.8))')         & 
		 (	                             & 
		      lgecmb(j)                  &
		    , tt      , Ppn  (j), tt, tt &
		    , Bpp  (j), Bpn  (j), tt, tt &
		    , Dpp  (j), Dpn  (j), tt, tt &
		    , dBpp (j), dBpn (j), tt, tt &
		    , dDpp (j), dDpn (j), tt, tt &
		    , d2Dpp(j), d2Dpn(j), tt, tt &
			, j=1,necmb                  &   
         )
	  close(1)
      print*, ' CMB coefficients are read!'
! ---------- end array for CMB  -------------

   end Subroutine ReadPCMB
!-----------------------------------------

!***************************************************
! Calculates all approximations to the             *
! coefficients  of e^+ + e^- and                   *
! nucleon photoproduction in CMBR  at z=0          *
!***************************************************
!  r = ln(E/eV)                                    *
!***************************************************
! c(1) = Bee, c( 2) = dBee, c( 3) = Ppn                 
! c(4) = Bpp, c( 5) = dBpp, c( 6) = Dpp, c( 7) = dDpp, c( 8) = d2Dpp  
! c(9) = Bpn, c(10) = dBpn, c(11) = Dpn, c(12) = dDpn, c(13) = d2Dpn  

   Subroutine AprBD(vlgE) 
   
	  use msimsl, only: DCSIEZ 
	  use comnoosc, only: r_0, ln10, FLIN 
	       
      implicit none

	  integer:: i

	  real(8), intent( in):: vlgE
	  
      real(8), dimension(1):: rv     ! lg(E/eV)
	  real(8), dimension(13,1):: res

	  rv(1) = vlgE


!**********************************************************************
      if(vlgE.le.lgee(1)) then
         res(1,1) = flin(vlgE,lgee(1),lgee(2),lgbetaee(1),lgbetaee(2))  ! lg B^ee
         res(2,1) = flin(vlgE,lgee(1),lgee(2),dbetaee (1),dbetaee (2))  ! dB^ee/dlnE 
      elseif(vlgE.gt.lgee(nee)) then
         res(1,1) = flin(vlgE,lgee(nee-1),lgee(nee),lgbetaee(nee-1),lgbetaee(nee))
         res(2,1) = flin(vlgE,lgee(nee-1),lgee(nee),dbetaee (nee-1),dbetaee (nee))
      else
         CALL DCSIEZ (nee, lgee, lgbetaee,1,rv,res(1,:) )
         CALL DCSIEZ (nee, lgee, dbetaee ,1,rv,res(2,:)	)
      endif
!**********************************************************************
      if(vlgE.le.lgecmb(1)) then
         res( 3,1) = flin(vlgE,lgecmb(1),lgecmb(2),Ppn  (1),Ppn  (2))
         res( 4,1) = flin(vlgE,lgecmb(1),lgecmb(2),Bpp  (1),Bpp  (2))
         res( 5,1) = flin(vlgE,lgecmb(1),lgecmb(2),dBpp (1),dBpp (2))
         res( 6,1) = flin(vlgE,lgecmb(1),lgecmb(2),Dpp  (1),Dpp  (2))
         res( 7,1) = flin(vlgE,lgecmb(1),lgecmb(2),dDpp (1),dDpp (2))
         res( 8,1) = flin(vlgE,lgecmb(1),lgecmb(2),d2Dpp(1),d2Dpp(2))
         res( 9,1) = flin(vlgE,lgecmb(1),lgecmb(2),Bpn  (1),Bpn  (2))
         res(10,1) = flin(vlgE,lgecmb(1),lgecmb(2),dBpn (1),dBpn (2))
         res(11,1) = flin(vlgE,lgecmb(1),lgecmb(2),Dpn  (1),Dpn  (2))
         res(12,1) = flin(vlgE,lgecmb(1),lgecmb(2),dDpn (1),dDpn (2))
         res(13,1) = flin(vlgE,lgecmb(1),lgecmb(2),d2Dpn(1),d2Dpn(2))
      elseif(vlgE.gt.lgecmb(necmb)) then
         res( 3,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Ppn  (necmb-1),Ppn  (necmb))
         res( 4,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Bpp  (necmb-1),Bpp  (necmb))
         res( 5,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dBpp (necmb-1),dBpp (necmb))
         res( 6,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Dpp  (necmb-1),Dpp  (necmb))
         res( 7,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dDpp (necmb-1),dDpp (necmb))
         res( 8,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),d2Dpp(necmb-1),d2Dpp(necmb))
         res( 9,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Bpn  (necmb-1),Bpn  (necmb))
         res(10,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dBpn (necmb-1),dBpn (necmb))
         res(11,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Dpn  (necmb-1),Dpn  (necmb))
         res(12,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dDpn (necmb-1),dDpn (necmb))
         res(13,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),d2Dpn(necmb-1),d2Dpn(necmb))
      else
         CALL DCSIEZ (necmb, lgecmb,   Ppn ,1, rv, res( 3,:) )
         CALL DCSIEZ (necmb, lgecmb,   Bpp ,1, rv, res( 4,:) )
         CALL DCSIEZ (necmb, lgecmb,  dBpp ,1, rv, res( 5,:) )
         CALL DCSIEZ (necmb, lgecmb,   Dpp ,1, rv, res( 6,:) )
         CALL DCSIEZ (necmb, lgecmb,  dDpp ,1, rv, res( 7,:) )
         CALL DCSIEZ (necmb, lgecmb, d2Dpp ,1, rv, res( 8,:) )
         CALL DCSIEZ (necmb, lgecmb,   Bpn ,1, rv, res( 9,:) )
         CALL DCSIEZ (necmb, lgecmb,  dBpn ,1, rv, res(10,:) )
         CALL DCSIEZ (necmb, lgecmb,   Dpn ,1, rv, res(11,:) )
         CALL DCSIEZ (necmb, lgecmb,  dDpn ,1, rv, res(12,:) )
         CALL DCSIEZ (necmb, lgecmb, d2Dpn ,1, rv, res(13,:) )
      endif

	  C(1) = 1d1**res(1,1)
      do i = 2, 13
	     C(i) = res(i,1)
	  enddo 

   end Subroutine AprBD
!**********************************************************************


!**********************************************************************
!     z1 = 1+z ;      r = ln(E/E_0)
!-----------------------------------------
! Cour(1) - U_r''  I 
! Cour(2) - U_r'   I  
! Cour(3) - U_r	   I 
! Cour(4) - Ppn_r  I 
! Cour(5) - Q_r	   I  
!-----------------------------------------
! c(1) = Bee, c( 2) = dBee, c( 3) = Ppn                 
! c(4) = Bpp, c( 5) = dBpp, c( 6) = Dpp, c( 7) = dDpp, c( 8) = d2Dpp  
! c(9) = Bpn, c(10) = dBpn, c(11) = Dpn, c(12) = dDpn, c(13) = d2Dpn  

   Subroutine AprCoefZ(z1,r)
	  use comnoosc, only: Lambda, Omega_m &
	                   , r_0, ln10       &
					   , gam             &
					   , evl              
	  implicit none

      integer:: i
 	  real(8), intent(in):: z1, r
	  real(8):: vlgE, z13, hz, fact_z, z1evl

	  z13    = z1*z1*z1
	  z1evl  = z1**evl 
	  hz     = dsqrt(Omega_m*z13 + Lambda)
	  fact_z = z13/hz

!++++++++++++++++++++++  Our World  +++++++++++++++++++++++++++
	  vlgE  = (r + r_0)/ln10 + dlog10(z1)   ! lg(E/eV) in our world	 

	  call AprBD(vlgE)
  
	  Cour(1) = 0.5d0*( C(6) + C(11) )
	  Cour(2) = 1d0 + C(1) + C(4) + 1.5d0*C( 6) + C( 7)            & 
                     + (C(9) + 1.5d0*C(11) + C(12) )  
	  Cour(3) = -2d0 + C( 1) + C( 2)                               &
	          + C(4) + C( 5) + C( 6) + 1.5d0*C( 7) + 0.5d0*C( 8)   &
	   + ( C(9) + C(10) + C(11) + 1.5d0*C(12) + 0.5d0*C(13) )
	  Cour(4) = 0d0
!--- Source generation function ---
! Q = (1+z)^3*(E/E_0)^{-gamma}
!	  if(r < 0) then
!         Cour(5) =	fact_z*dexp(-2d0*r)
!	  else
	     Cour(5) =	fact_z*dexp(-gam*r)*z1evl
!	  endif
!----------------------------------
      do i = 1,4
         Cour(i) = fact_z*Cour(i)
	  enddo


   end Subroutine AprCoefZ
!**********************************************************************

end module coeffnoosc
